Phase separation of a model binary polymer solution in an external field. 
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The phase separation of a simple binary mixture of incompatible linear polymers in solution 
is investigated using an extension of the sedimentation equilibrium method, whereby the osmotic 
pressure of the mixture is extracted from the density profiles of the inhomogeneous mixture in a 
gravitational field. In Monte-Carlo simulations the field can be tuned to induce significant inho- 
mogeneity, while keeping the density profiles sufficiently smooth for the macroscopic condition of 
hydrostatic equilibrium to remain applicable. The method is applied here for a simplified model of 
ideal, but mutually avoiding polymers, which readily phase separate at relatively low densities. The 
Monte-Carlo data are interpreted with the help of an approximate bulk phase diagram calculated 
from a simple, second order virial coefficient theory. By deriving effective potentials between poly- 
mer centres of mass, the binary mixture of polymers is coarse-grained to a "soft colloid" picture 
reminiscent of the Widom-Rowlinson model for incompatible atomic mixtures. This approach sig- 
nificantly speeds up the simulations, and accurately reproduces the behaviour of the full monomer 
resolved model. 

PACS numbers: 



I. INTRODUCTION 

Polymer solutions segregate into polymer-rich and 
polymer-poor phases when the temperature is lowered 
below the theta-point. Polymer blends nearly always 
demix in the melt, because the connectivity constraints 
strongly reduce the entropy of mixing 1 . Less is known, 
however, about the phase behaviour of dilute or semi- 
dilute polymer solutions involving two or more poly- 
meric species. In particular, the quality of the solvent 
is generally not the same for different polymeric species, 
which may lead to competing mechanisms in the sub- 
tle free-energy balance between entropic and energetic 
contributions. In the case of binary mixtures of incom- 
patible polymers in a common good solvent, excluded 
volume effects lead to strong deviations from mean-field 
behaviour 2,3 . In this paper we consider a highly sim- 
plified model of mutually repelling polymer in a com- 
mon solvent. The system under consideration extends 
the familiar Widom-Rowlinson model^ for incompatible 
atomic mixtures to polymer blends in solution. Not sur- 
prisingly, the mutual incompatibility of the two species 
in the athermal model leads to demixing above a critical 
density. The additional feature introduced in the present 
work is the presence of an externally applied field, specif- 
ically a gravitational field, acting on the polymers, which 
induces a spatial inhomogencity. 

The interest of the applied field is two-fold. The first is 
essentially methodological: as we have shown recently for 
cases of single-component polymer solutions^ and of di- 
block copolymer solutions^, the density profiles charac- 
terising the inhomogeneous distribution of polymers may 
be exploited to extract the osmotic equation-of-state over 
a wide range of concentrations in a single simulation*. A 
second motivation is that applying an external field al- 
lows us to rapidly search through phase space to locate 
phase-transitions, although it must be kept in mind that 



adding a field affects the phase behaviour, as demon- 
strated in a recent theoretical investigations of colloid- 
polymer mixtures^. 

Finally, we demonstrate in this paper the merits of 
coarse-grained representations of complex polymeric sys- 
tems, whereby the total number of degrees of free- 
dom is drastically reduced by tracing out individual 
monomer coordinates to determine effective interactions 
acting between polymer centres of mass. This strat- 
egy has proved very successful for solutions of linear 
polymers2iifiiii, mixtures of asymmetric star polymers**, 
di-block copolymers^ and is extended here to linear poly- 
mer mixtures, in an effort to map out their phase dia- 
gram. Our coarse-graining strategy is similar in spirit to 
that successfully applied by the Klein group to surfac- 
tants and other self-assembling systems^. 

II. MODEL AND METHODOLOGY 

Consider a binary polymer mixture of linear homopoly- 
mers with Ma and Mb monomers respectively. For com- 
putational purposes the two species are assumed to " live" 
on a cubic lattice of L x x L y x L z sites. Monomers oc- 
cupy lattice sites, and the lattice spacing coincides with 
the segment length b. Let Na and Nb be the total num- 
bers of polymers of the two species. The overall poly- 
mer densities (coils per unit volume) are p a = N a /V, 
where V — L x L y L z b 3 (we set b to 1 in this paper). 
An external force field along the z-direction, acting on 
monomers of each species a, and deriving from the po- 
tential ip a (z)(a = A or B) will induce an inhomogene- 
ity in the mixture characterised by two density profiles 
p a (z). Alternatively, one may define the overall polymer 
density p(z) — pa{z) + Pb(z) and the concentration ratio 
profile: 

x{z) = p A (z)/p(z) (1) 
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The density profiles are normalised such that: 



p a (z)dz = n a 



L X Ly 



(2) 



If the inhomogeneity induced by the external field varies 
sufficiently slowly in space, the binary mixture obeys the 
hydrostatic equilibrium condition: 



dP(z) , d^ A (z) ,d^ B {z) 
= -pA(z)— — Pb{z)- 



dz 



dz 



(3) 



where P(z) is the local osmotic pressure of the binary 
polymer mixture in solution. The macroscopic condition 3 
may be shown to follow from the "local density approxi- 
mation" (LDA) within density functional theory of inho- 
mogeneous fluids^. If the profiles pa(z) an d Pb( z ) are 
computed in Montc-Carlo (MC) simulations of binary 
mixtures subjected to the external field, P(z) may be 
determined by integration of equation [3] The equation- 
of-state of the inhomogeneous mixture along an isotherm, 
P(pA, Pb) or P(p, x) may then be extracted by eliminat- 
ing the altitude z between P(z), pa(z) and pb(z) (or p(z) 
and x(z)). 

In practice, the external field is chosen to be the grav- 
itational field acting on the centre-of-mass of the poly- 
mers: 



ip a (z) = M a m a gz 



(4) 



where m a is the buoyant mass of monomers of species 
a and g is the acceleration due to gravity. Since gravity 
acts as a confining field, it is convenient to take L z — » oo. 
Although, under normal experimental conditions, gravity 
is insufficient to induce significant inhomogeneity in a 
polymer solution, (ultra-centrifugation would be required 
to induce a measurable effect) here gravity may be tuned 
at will in simulations to achieve sedimentation lengths 
C = kBT/m a g comparable to radii of gyration R ga , and 
hence a significant variation of the density profiles. 

Substitution of equation0]into equation[3]and integra- 
tion leads to : 

/>oc />oo 

P(z)/k B T = 1/Ca J Pa (z')dz' + l/C/3 J pp{z')dz'.{5) 

The equation-of-state P(p a ,pp) of the heterogeneous 
mixture may thus be determined from the density pro- 
files Pa(z) and pp(z) as explained above. In practice MC 
runs are carried out for several overall concentrations, in 
order to map out the complete equation-of-state. It is 
important to realise that, in general, not only the overall 
density p(z), but also the concentration x{z) varies with 
altitude, as may be easily verified for the simple case of 
a binary mixture of ideal (non-interacting) polymers. 

We have applied the hydrostatic equilibrium method 
to a very simple non-trivial model of a binary polymer 
mixture which demixes in solution. In the model A and 
B polymers freely penetrate coils of the same species, 
i.e. follow random walk statistics, while A-B pairs ex- 
perience excluded volume between their monomers i.e. 



behave like mutually avoiding walks. The model is one 
of a family of six XYZ models, where X, Y, and Z la- 
bel the statistics of A-A, A-B and B-B interaction pairs. 
Thus each of X, Y, and Z can take one of two "values" , 
I for ideal (non-interacting) and S for self (or mutually) 
avoiding. In particular SSS refers to a binary mixture of 
polymers in which all interactions are self avoiding, i.e. a 
homogeneous solution of SAW polymers, while the model 
under consideration here is the ISI model. Note that 
the ISI model generalises the familiar Widom-Rowlinson 
model^ to polymer mixtures; a related block co-polymer 
model, where A and B are tethered, has recently been 
investigated by our group 6 . While the latter leads to 
micro-phase separation, the present binary ISI model is 
expected to undergo bulk phase separation due to the in- 
compatibility of the A and B components. For the sake 
of simplicity, the following investigations are restricted 
to the case of equal length polymers (Ma = Mb — M). 
Moreover, without loss of generality, the monomer masses 
may be assumed to be equal (itia = tub = m), so 
that (a = Cb — C Note that this model is athermal. 
Polymers were subjected to pivot and translational MC 
moves, typically 1 x 10 6 — 1 x 10 s MC moves per polymer. 



III. COARSE-GRAINED DESCRIPTION 

Following a known strategjaiii2*ii, one can calculate 
effective interactions v a p(r) between CMs of A-A, A-B 
and B-B pairs, by averaging over the conformations of 
two polymers for fixed values of the distance r between 
their CMs. Explicitly, for an isolated pair of A and B 
polymers: 



v a p(r) = -k B TlnP A B(r) 



(6) 



where Pab{t) is the probability (averaged over all ran- 
dom walk configurations of the two polymers) that there 
is no overlap between monomers of A and monomers of 
B when their CMs are held at a distance r. Obviously 
VAA{f) = vbb(t) = for the ISI model, since like-species 
overlaps are allowed. Effective pair potentials like vab 
in eq. are easily computed in MC simulation, as ex- 
plained on more detail in the referencesiSiii. Examples 
are shown in figure ^ for the ISI, ISS and SSS mod- 
els; since the two species are of equal length, the SSS 
model is equivalent to a one-component system of self- 
avoiding walk (SAW) polymers, which has been thor- 
oughly investigatedi2iii*ii While the potentials vab{t) 
for the SSS and ISS models are very close, the effective 
potential for the ISI model is considerably more repul- 
sive,with a value at full overlap, vab(t — Q) — 3.45/csT. 
This enhancement of the effective repulsion may be eas- 
ily understood by noting that the A and B coils are more 
compact (their radii of gyration are those of ideal chains, 
which scale as R g ~ M 5 , as compared to chains with 
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SAW internal statistics, which scale as R g 
Hence the overlap probability is higher (i.e. Pab{t) is 
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FIG. 1: Effective potential vab(t), in units of fcgT, between 
the centres of mass of two polymers for ISI, ISS and SSS 
models. These potentials were generated at zero density using 
Monte Carlo simulations of two L = 500 polymers on a simple 
cubic lattice. 



lower) in the former compared to the latter case, for any 
given r/R e ff where the effective radius is defined by 

R 2 e }f = \{R g A 2 + R gB 2 ). (7) 

The second virial coefficient for the ISI model may be 
calculated from Vab{t) according to: 
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B AB = -27T / ( e -^B(r)/fe B T _ 1)r 2 dr (8) 
JO 

and turns out to be Bab — 9.5i?^ for our L = 500 ISI 
model lattice polymers. 

It has been showniS that the effective pair poten- 
tial between linear SAW polymers varies significantly as 
the polymer density is increased from the infinite dilu- 
tion limit, particularly so beyond the overlap density 
p* = 3/4:TrRg. This effect is expected to be less im- 
portant for the ISI model, since it will be shown that 
phase separation, which leads to a dramatic reduction of 
the number of A-B contacts, sets in for p greater than 

~ IP ■ 

IV. EQUATION-OF-STATE OF THE 
SYMMETRIC MIXTURE 

In order to determine the full equation-of-state 
P(Pa,Pb) or P(p,x) by the hydrostatic equilibrium 
method presented in section 2, simulation of the ISI 
model must be carried out on systems involving sev- 
eral overall compositions. In the equimolar case, the 
ISI model is fully symmetric with respect to the A and 
B species, so that the local concentration x(z) remains 
constant (equal to x = |) at all altitudes. Examples of 



FIG. 2: The centre of mass and monomer density profiles for 
the symmetric equimolar ISI model, with the same profiles 
on a logarithmic scale shown as an insert. These profiles were 
generated from simulations of 1600 polymers of length L = 
500 in a box with a base of L x = L y — 200 that is open 
in the positive z direction, along which the density profiles 
are plotted. The dashed line on the logarithmic scale (inset) 
shows the exponential behaviour of ideal polymers p(z) ~ 
exp(— a/C), with the gravitational length £ = 2.19_R S . 



MC-generated monomer and CM density profiles pa(z) 
and pb{z) of the equimolar ISI model are shown in figure 
El As expected, the density profiles of the two species 
are identical within statistical errors. The two profiles 
exhibit a depletion zone (which is more pronounced for 
the CM profile) near the bottom, followed by a peak at 
z ~ 2-R ff , which is again more pronounced for the CM 
profile. At higher altitudes the profiles go over to the 
exponential barometric behaviour of ideal particles, as 
clearly seen in the inset. The intermediate behaviour 
signals phase separation (as discussed below). Beyond 
the peak, the CM and monomer profiles are virtually in- 
distinguishable. Only the latter part is used in the inver- 
sion procedure to extract the osmotic equation-of-state, 
since the rapid variation at low altitudes is incompatible 
with the LDA, and hence with the macroscopic condi- 
tion of hydrostatic equilibrium embodied in eq. [S] The 
resulting equation-of-state Z = P/pksT is plotted in fig- 
ure After an initial linear increase, Z goes through 
a maximum around p/p* < 0.5, before decreasing and 
apparently saturating at a value above one at high den- 
sities. This unusual behaviour may be traced back to to 
the phase separation, which is evident in the snapshot 
of a typical configuration shown in fig 0] and which will 
be analysed in more detail in section 5. The existence of 
phase boundaries leads one to expect a substantial de- 
pendence of the measured equation-of-state on the size 
of the simulation box. 

To illustrate this size dependence, we have carried out 
MC simulations of the coarse-grained representation of 
the ISI model, based on the effective potential vab(?) m ~ 
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FIG. 3: The equation-of-state Z = PP/ p for different box 
sizes, denoted by R g /LxY, calculated by coarse grained simu- 
lations and compared, at R 9 /Lxy ~ 0.0365, to an equivalent 
fully resolved monomer-level simulation 



The small discrepancies are probably due to the fact that 
we have neglected the density dependence of the effective 
potential vab{t), for which we have used the zero den- 
sity result shown in figure ^ The simulation data for the 
coarse-grained model plotted in figure [3] clearly illustrate 
the size dependence at high densities, i.e. in the region 
of phase coexistence. This is expected since the smaller 
the system, the larger the interfacial contribution to the 
equation-of-state; the former will become negligible only 
in the thermodynamic limit (R g / \L xy) ~^ 0). 



V. INTERPRETATION OF THE SIMULATION 
DATA 

The equation-of-state data of figure [3] can be under- 
stood in terms of an underlying phase separation, as sug- 
gested by fig 0] The phase diagram of the ISI model 
can be calculated, at least approximately, from the virial 
expansion of the free energy. The free energy per unit 
volume can be expanded as: 




FIG. 4: Snapshot of a configuration for Na = Nb = 800 
polymers, showing the CM of each as a sphere of radius R g . 
The sedimentation lengths are (a = Cs = 2.11_R 9 , and the 
open box has a base of length Lx = Ly = 200. Note that 
the existence of two interfaces is due to the use of periodic 
boundary conditions in the x and y directions. 



troduced in section 3 (fig^) for three values of the ratio 
Rg/(LxY) (where Lxy = Lx = Ly is the length of the 
square base of the simulation cell). These were carried 
out for Na = Nb = 800 — 2000 coarse-grained particles, 
with a sedimentation length C, — 2R g . Such simulations 
are typically two orders of magnitude faster than those 
based on the full monomer-level ISI model for the same 
number of polymers. For R g / (Lxy) = 0.0365, the agree- 
ment between the equations of state calculated with the 
two representations is seen to be satisfactory in figure 
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ft + Ib + 2ar(l ~ x)B ABP 2 + 0(p 3 ) (9) 



where = p a (ln(N^p a ) — 1) and Bab is the second 
virial coefficient of equation (J8J. Differentiation of @ 
with respect to p yields: 



Z = 



P 



pk B T 



= df/dp -f/p=l + 2x(l - x)B AB p + O(p 2 \10) 



Above a critical density p c , the free energy @ leads to 
phase separation; by symmetry the critical concentra- 
tion c c = 1/2, while p c = I/Bab- Inserting the nu- 
merical value Bab = 9.5Rg reported in section 3, one 
finds Pel P* = 0.44. Following an analysis similar to one 
presented in ref^ one can calculate the spinodal curve 
analytically: 
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(11) 



The binodal is easily calculated numerically, and results 
are summarised in figure[3 segregation of the two species 
is seen to be nearly complete for p/p* ~ 1. The equation- 
of-state data of figure|21may now be "read" as follows. As 
the altitude z decreases in the polymer sediment, the lo- 
cal density p(z) increases (c.f. figureEJ. Up to p/p* < 0.4 
the system remains in the one phase region of the phase 
diagram in figure[Sl and the equation-of-state Z increases 
linearity with p/p* according to (|10|l : the slope of the lin- 
ear portion of the MC-generated equation-of-state agrees 
well with that predicted by JUJJ (with Bab — 9.5R g ). 
When the local density p(z) increases beyond the criti- 
cal density p c /p* — 0.44, the system enters the phase- 
coexistence region bounded by the binodal curve and 
separates symmetrically into two mixtures with compo- 
sitions x and 1 — x which may be read from figure |S] 
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FIG. 5: Spinodal and binodal demixing lines for the ISI 
model, calculated from second virial coefficient expansion of 
the free energy described in Eqs. 1^1- 1111 
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FIG. 6: The equation-of-state from the largest coarse-grained 
system simulations is compared to theoretical curves, calcu- 
lated from the binodal and the spinodal lines of a second virial 
expansion of the free energy described in Eas.l^l- llll 



The excess contribution to the pressure of the coexisting 
phases decreases with density, because the number of A- 
B contacts drops sharply in the A and B-rich phases as 
the overall density rises. At densities p > p*, the nearly 
complete segregation means that coexisting phases be- 
have practically as ideal cases of A or B polymers. Hence, 
for sufficiently large systems such that the interfacial con- 
tribution to the local stress is negligible compared to the 
bulk contributions, the equation-of-state Z must return 
to its ideal gas value 1 for p ^> p* . 

The theoretical equation-of-state is compared in figure 
to the MC data for the largest system. The agreement 
is seen to be only qualitative. The theoretical curve has 
a cusp point (rather than a rounded maximum) at the 



density where phase separation starts. In the simula- 
tions, the system may enter the metastable region com- 
prised between the binodal and spinodal curves in figure 
[S] In that case the mixture may reach higher densities be- 
fore segregating. The equation-of-state calculated on the 
assumption that the mixture separates into two phases 
along the spinodal curve is also shown in figure [5] It is 
gratifying to realise that the simulation data lie between 
the binodal and spinodal equations-of-state, which con- 
stitute two extreme limits of the actual scenario. 
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FIG. 7: Density and concentration profiles Pa{z), pb(z), x(z) 
(inset) of a 70:30 mixture of ISI polymers. 

In the more general case where the overall composition 
of the binary polymer solution x ^ 0.5 (i.e. Na ^ Nb), 
the local concentration ratio, x(z) will no longer be con- 
stant, and the extraction of Z(p, x) from equation JSJ re- 
quires several simulations of systems with different com- 
positions. No attempt is made at a full investigation in 
the present paper, but the situation is illustrated in fig- 
ure [7| where density and concentration profiles obtained 
from MC simulations of the full monomer description of 
a system with Na = 480 and Nb = 1120 (a; = 0.3) poly- 
mers of La — Lb = 500 in an open square based box 
of side length Lxy — 200 are plotted. The concentra- 
tion profile x(z) is seen to vary substantially at altitudes 
where phase separation is expected to take place. 

In this case, as well as for lower overall concentration 
fractions of A particles, we observe the formation of a ver- 
tical, cylindrical interface, rather than two planar inter- 
faces as in figure 0] This observation can be understood 
as the consequence of a minimisation of the area over 
which unfavourable A-B contacts take place. Assuming 
that the distributions pa(z) and pb{z) of the two species 
are identical irrespective of the geometry of the interface, 
and that the height H of the two types of interfaces is the 
same, and noting that the radius r of the cylinder is re- 
lated to the concentration x by x = n(r / Lxy) 2 , one finds 
that the areas of the two planar and cylindrical interfaces 
are identical when x = l/V; below this value, the cylin- 
drical interface will be preferred. We explicitly checked 
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this by performing simulations for Na + Nb = 1600 poly- 
mers at x = 0.1, 0.2, 0.3, 0.4 and 0.5. As expected, for the 
lowest three concentration ratios we find cylindrical in- 
terfaces, and for the two higher concentrations we find 
planar interfaces. 

VI. CONCLUSION 

We have extended the sedimentation equilibrium 
method for the computation of the osmotic equation-of- 
state of polymer solutions to the case of binary poly- 
mer mixtures undergoing phase separation. The method 
was applied to the athermal ISI model, where polymer 
species A and B (assumed here to be of the same length) 
behave like ideal (random walk) polymers, but are mutu- 
ally avoiding, i.e. monomers of different species cannot 
occupy the same lattice sites. This incompatibility leads 
to segregation for overall polymer densities higher than 
about one half of the overlap density p* . This segrega- 
tion leads to a somewhat unusual variation of the osmotic 
pressure with density, since outside the interfacial region, 
the mixture behaves as an ideal solution both in the low 
and high density limits. 

The MC generated equation-of-state agrees well with 
the data obtained from a coarse-grained representation 



based on an effective interaction between the CMs of A- 
B pairs, and agrees semi-quantitatively with the predic- 
tions of a simple analytic theory based on a second virial 
expansion. The cusp in the equation-of-state predicted 
by the latter is rounded into a maximum in the MC data 
due to finite size effects. The present investigation sug- 
gests that the phase behaviour of segregating mixtures 
can be extracted from a measurement of the density pro- 
file of the inhomogeneous mixture in an external field. 
Although our model is highly oversimplified, we have 
demonstrated that our coarse-graining strategy can be 
applied in a straight-forward fashion to the demixing of 
polymer solutions. We plan to extend this methodology 
to more realistic polymer models, and to grand canonical 
ensemble simulations. 
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